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Abstract 

This study is motivated by problems related to environmental transport on river 
networks. We establish statistical properties of a flow along a directed branching 
network and suggest its compact parameterization. The downstream network trans- 
port is treated as a particular case of nearest-neighbor hierarchical aggregation with 
respect to the metric induced by the branching structure of the river network. We 
describe the static geometric structure of a drainage network by a tree, referred to 
as the static tree, and introduce an associated dynamic tree that describes the trans- 
port along the static tree. It is well known that the static branching structure of 
river networks can be described by self-similar trees (SSTs); we demonstrate that 
the corresponding dynamic trees are also self-similar. We report an unexpected 
phase transition in the dynamics of three river networks, one from California and 
two from Italy, demonstrate the universal features of this transition, and seek to 
interpret it in hydrological terms. 

1 Introduction 

The topology of river networks has been extensively studied over the past decades and 
stream ordering schemes, as well as statistical self-similarity concepts, have been explored 
to a considerable extent [see Norton, 1945; Strahler, 1957; Shreve, 1966; Tokunaga, 1978; 
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Mandelbrot, 1983; La Barbera and Rosso, 1987; Marani et al, 1991; Rodriguez-Iturbe et 
al, 1992; Peckham, 1995; Badii and Politi, 1997; Rodriguez-Iturbe and Rinaldo, 1997; 
Turcotte, 1997; Sposito, 1998; Peckham and Gupta, 1999; Pelletier and Turcotte, 2000; 
5uro? a/., 2000; Dodds and Rothman, 2000; da Costa et al, 2002; and references therein] 

What has been less studied, however, is how the static topology of a river network 
affects and is affected by the dynamical processes operating over this network. For 
example, consider a directed tree that represents a river network, and assume we are 
interested in the mixing of water, solutes and sediments as they move downstream in 
reaches of variable lengths and merge at junctions of the river network. One might want 
then to attach a "metric" to the nodes of this tree, such as the distance to the nearest 
source, and consider the notion of a dynamic tree, superimposed on the template of the 
underlying static tree. 

This dynamic tree is likely to have a different hierarchy and topology than the static 
one. Depending on the dynamics, for example, some of the static-tree branches might be 
completely cut off, either due to a blockage that prevents transport along these branches 
or due to the absence of conditions to generate sediment or nutrient for downstream 
transport. In this case, the dynamic tree will have a different hierarchy than the static 
one, and this difference might affect the scaling of fluxes that participate in defining the 
envirodynamics on the network of interest. In general, a static tree of a given florton- 
Strahlcr order could become a dynamic tree of a lesser or higher order, depending on the 
superimposed dynamics. 

The purpose of this paper is to study the dynamic topology of directed trees, starting 
with several simple cases, first synthetic and then realistic. The direction we use is obvi- 
ously "downstream," i.e., from the leaves to the root of the tree. We focus on a dynamic 
hierarchy built on the concept of "connectivity" : once two streams are connected, they 
both influence the downstream dynamics. One can thus imagine that two order- 1 streams 
of different lengths, li and I2, merge at a node but do not automatically give rise to an 
order-2 stream, as would be the case in the standard Horton-Strahler ordering scheme; 
instead, we keep track of length and the assigned order becomes 2 only when the running 
index of length becomes max(Zi, I2). 

Alternatively, one might keep track of time, rather than length: the two are equivalent 
if the flow velocity is constant along all the branches, which we will assume in the present 
paper, for simplicity's sake. In other words, a dynamical node of order 2 is created 
only when the fluxes from both order- 1 streams do reach the connecting node. Such 
considerations will result in a different ordering of the dynamic tree compared to the 



2 



static one. Moreover, the newly created dynamic tree will be time- oriented, a property 
that is absent in conventional static trees. 

We approach the problem of hierarchical djTiamics of river networks using general 
concepts of hierarchical aggregation, which studies how multiple individual particles 
(molecules, species, individuals, etc.) merge (aggregate, collide) with each other to form 
clusters in different physical, chemical, biological, or sociological settings. A major role 
in such studies is played by the notion of cluster dynamics. This concept refers to the 
situation when a system that contains an infinite number of interacting particles can be 
decomposed into finite clusters that move independently of each other for some random 
interval of time. After this time, the particle interactions give rise to infinite-range cor- 
relations, although the system can be decomposed into another set of finite independent 
clusters, and so on. 

In the 1970s, Ya. G. Sinai developed a self-consistent mathematical formalism and 
proved the existence of cluster dynamics for some particle systems in statistical mechanics 
[Sinai, 1973, 1974]. The ideas of cluster dynamics have been applied to plasma physics, 
economics, and the study of precursory patterns for extreme events in geophysics [Rotwain 
et ai, 1997; Molchan et ai, 1990; Keilis-Borok and Soloviev, 2003]. Recently, Gabrielov 
et al. [2008] evaluated numerically the cluster dynamics of elastic billiards, leading to the 
detection of what appear to be the first genuine phase transitions and scaling phenomena 
that develop in time, rather than with respect to a control parameter, such as temperature 
T or density; i.e., a transition occurs and scaling develops as time t evolves toward a 
critical value t* , rather than as the parameter T crosses a critical value T*. 

In this paper, we adapt the concept of cluster dynamics to environmental transport 
on river networks. Notably, we obtain a remarkably similar, and equally unexpected, 
phase transition in the cluster dynamics of river networks and attempt to interpret it in 
this context. We also study the statistical properties of the dynamic trees introduced 
herein. It is well known that the static branching structure of river networks can be 
described by self-similar trees (SSTs); we demonstrate, using three actual river basins, 
that the corresponding dynamic trees are also self-similar. 

This paper is structured as follows. We review in Section [2]the terms and concepts rel- 
evant to the hierarchical analysis of branching structures, including the Horton-Strahler 
and Tokunaga branching taxonomies. Section [3] introduces the concept of a dynamic tree 
that is associated with a given static tree, by using two examples from river transport. 
Section m describes two types of static trees analyzed in this study. The first type reflects 
the well-formed "river network" of a basin. The second type reflects the "unchannelized 
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drainage network" ; this network is composed of drainage paths that are not permanent 
channels but are perpendicular to the topographic contour lines and follow the steep- 
est downstream gradient. In other words, this drainage network is formed by paths of 
"zero-order" basins or hillslopes. Hierarchical aggregation is described in greater depth, 
and with additional examples from several fields, in Section |U along with an abstract 
metric space setup. Three actual river networks, from California and Italy, are analyzed 
in Section [6l A summary and discussion follow in Section [7l 

2 Main concepts and definitions 

This section introduces the main concepts used in the analysis of branching structures, 
along with their definitions and illustrative examples. 

2.1 Trees 

A tree T is a set of nodes connected by vertices (also called edges or links) in such a 
way that there are no loops, i.e. there are no closed paths formed by distinct edges 
(see Fig. [1]). A rooted tree has one special node designated as a root. In a rooted tree 
each connected pair of nodes has a parent-child relationship, with the parent being the 
element that is closest to the root [Athreya and Ney, 1972]. The nodes with no children 
are called leaves. The depth di of a node z in a rooted tree is defined as the number of 
edges between this node and the root. The depth D of a tree is the maximum of the 
depths di over all the leaves /. 

In this study we will work with binary trees. In a binary rooted tree, each node may 
have either two or no children. This means that each internal node i (every node except 
for the root and the leaves) is connected to three other nodes: one is a parent of and 
the other two are its children. The root is only connected to two children, and each leaf 
is connected to a single parent. A complete binary tree is a rooted tree such that all its 
leaves have the same depth. Our interest for binary trees is motivated by the observation 
that many natural phenomena exhibit binary branching. For example, in river networks, 
it is unlikely for three or more streams to merge at exactly the same point, while in gas 
dynamics it is unlikely that more than two molecules will collide at the same time. 

In our study of river transport, the tree T will represent a drainage network; see 
Fig. [H Hence the nodes correspond to the merging points of streams and vertices to the 
stream segments between these points, while the network's sources are the leaves, and 
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the outlet is the root of the tree. 



2.2 Branching-order taxonomies 

In many apphcations, there is a need to order the nodes according to their importance 
in forming the entire hierarchy; this importance often corresponds also to relative size. 
For instance, in a botanical tree the leaves are the most delicate, smallest elements; the 
intermediate levels are formed by consecutively wider branches, while the most heavy, 
robust element of the plant is its trunk. Likewise, one naturally distinguishes in a river 
network between minor and major tributaries, according to the amount of water that 
they are able to carry. 

In a complete tree, the node ordering task is quite straightforward since a node's 
order can be chosen to be inversely proportional to its depth: "the deeper, the smaller" . 
The problem, however, becomes more complicated when one deals with an incomplete 
tree; in this case, the depth can no longer serve as a proxy for size, since the leaves, while 
being the smallest elements, will often be assigned indices that are as large as those of 
much heavier internal nodes. 

Horton [1945] developed a convenient way to order hierarchically organized river trib- 
utaries; this method was later refined by Strahler [1957] and further expanded by Toku- 
naga [1978]. Currently, the so-called Horton-Strahler and Tokunaga ordering schemes 
are standard tools of branching analysis. 

2.2.1 Horton-Strahler ordering 

Each leaf in a binary rooted tree is assigned a Horton-Strahler (HS) order r(leaf) = 1; 
see Fig. [2^. Each node p, which is the parent of nodes Ci and C2, is assigned a Horton- 
Strahler order r{p) according to the following rule [Horton, 1945; Strahler, 1957; Newman 
et al, 1997]: 



A branch is defined as a union of connected nodes with the same order. We will 
denote by A^,. the total number of branches of order r. Notice that each branch has linear 
structure: two children of the same parent can not belong to the same branch. 

In a tree with n leaves, the longest branch can be formed by [n — 1) nodes; this is 
the case when two leaves merge together to form an order-2 branch and then all other 
leaves join this branch one by one. We refer to this situation as exhaustive branching. It 





5 



is readily seen that each leaf is always an order- 1 branch. An order-2 branch is created 
by merging two leaves and can consist of more that one node, depending on the leaves 
that join it; an order-2 branch that consists of two nodes is highlighted in Fig. [2^. The 
order Q of a. tree is the maximal order of its branches (or nodes). 

In a complete tree, each branch consists of a single node since the children of an 
order-r node always have the same order (r — 1). In such a tree, the HS order is uniquely 
determined by the node depth d via r = D — d + 1, where the tree depth is D = Q. 

2.2.2 Tokunaga indexing 

Tokunaga indexing [Tokunaga, 1978; Peckham, 1995; Newman et ai, 1997] extends 
upon the Horton-Strahler orders; it is illustrated in Fig. [2]d. This indexing focuses on 
incomplete trees by cataloging the merging points between branches of different order. 
A first-order branch that merges with a second-order branch is indexed by "12" and the 
total number of such branches is denoted by A first-order branch that merges with 
a third-order branch is indexed by "13" and the total number of such branches is N13, 
and so on. In general, Nij for j > i denotes the total number of order-z branches that 
join an order-j branch. 

The Tokunaga index Tij is the number of branches of order i that merge with a branch 
of order j, normalized by the total number of branches of order j; in other words, is 
the average number of branches of order i < j per branch of order j: 



T- • = ^ (2) 



Merging of branches of different orders is referred to as side branching. It is easily 
seen that side branching is absent in a complete tree, and "a tree with side branching" is 
synonymous to "an incomplete tree." For incomplete trees, the side-branching indices be- 
come increasingly important as they help to define a tree's structure, possibly indicating 
properties that are unique to specific classes of trees. 

For consistency, we denote the total number of order-z branches that merge with other 
order-i branches by Nu and notice that in a complete binary tree Na = 2A^j+i. This 
allows us to formally introduce the additional Tokunaga indices: 

T., = ^ = 2. 

The set {Tij : I < i, j < fl} of Tokunaga indices provides a complete statistical descrip- 
tion of the branching structure of an order-f2 tree. 
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2.2.3 Other node statistics 

We introduce here two node statistics relevant to our river transport study: 

• the number of nodes (or links) within a branch i is denoted by cf, and 

• the magnitude rrii of a node i is the number of leaves that descend from i; in other 
words, the magnitude of a branch is the number of the sources upstream of it. 

Magnitude measures the complexity of the river structure upstream from a given 
branch. We notice that each leaf (source) has unit magnitude, mieaf = 1, and the 
magnitude of a parental node p is the sum of the magnitudes of its children ci and C2'- 

nip = rric-, + rric^. (3) 

Accordingly, a node or order r has magnitude m > 2^^^, with equality being attained 
only for a complete binary tree. The average number of nodes and average magnitude of 
an order-r branch are denoted by Cr and Mr respectively. 



2.3 Self-similar trees 

The concept of self-similarity provides a powerful tool for describing and studying trees. 
A self- similar tree (SST) is defined by the constraint 

Ti^i+k = Tk for k = 1,2,.... (4) 

E. Tokunaga was probably the first to study SSTs, and considered an additional 
constraint on the branching indices [Tokunaga, 1978]: 

%ti = c, or Tk = a c^'^ for a, c> 0. (5) 

The SSTs that satisfy dSj) are called Tokunaga trees. 



2.4 Horton laws 

Empirically, the average values of branching statistics for the observed river basins depend 
exponentially on the order r: 

Nr = NoR'l-\ (6) 
Mr = R'm\ (7) 

Cr = Co Rq (8) 
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for some positive constants A^o and Cq. Such relationships are called Horton laws; the 
bases Rb,Rm, and Rc of the exponential relatonships are called stream ratios. 

McConnell and Gupta [2008] showed that the Horton laws ([6]), (JTj) hold asymptotically, 
i.e. for r oo, in a self-similar Tokunaga tree; they also proved that Rb = Rm- 
Moreover, Zaliapin [2009] demonstrated the stream ratio inequality 

Rb = Rm < Rc, (9) 

that had been conjectured by Peckham [1995]. In addition, Zaliapin [2009] demonstrated 
that the Horton laws hold, under some additional assumptions on the Tokunaga indices 
Tfc, for self-similar trees that do not necessarily satisfy the Tokunaga condition ([5]). 

3 Static vs. dynamic trees: Network envirodynam- 
ics 

The topological structure of a river network is well described by a tree, which we denote 
by Ts and call the static tree. To describe the downstream transport on Ts we now 
introduce a dynamic tree T/j, which can be interpreted as follows. Imagine that we 
inject a dye simultaneously into all the sources of our river network, represented by the 
leaves of T5, and the dye starts propagating down the river, from the sources to the 
outlet, with the same constant velocity along all the streams. The tree Td describes the 
time-dependent history of the mergings of the colored streams. 

Next, we consider two detailed examples that will clarify this important concept. We 
restrict ourselves to the simplest case of constant velocity along all the streams; taking 
this velocity to be unity, time and length scales can be interchanged. An extension 
to spatially or temporally variable velocities is straightforward: we shall see that the 
dynamic tree Td is completely determined by the static tree Ts and the set of time 
delays Ti necessary for the dye to propagate from a node i to its parent. 

3.1 Synthetic example 

Figure [3] shows how to construct the dynamic tree for a basin with four sources a, b, c, 
and d. The static tree for this basin is a complete binary tree shown in the top right 
panel. The same tree with the link lengths explicitly shown is placed in the top row of 
panels; the top left panel indicates the values of these lengths. 
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The consecutive phases of construction of the dynamic tree are shown in the bottom 
row of panels. At step (the leftmost top and bottom panels), all the links in the tree 
are "empty" (dashed lines) and the dye is injected into the sources a, b, c, and d. 
Accordingly, we have four disconnected clusters of colored flux; they correspond to four 
disconnected nodes in the lower left panel. Each step in the figure is a snapshot of the 
process after a unit time interval; recall that we only use constant velocity in this paper 
and, without loss of generality, this velocity equals unity. 

At step 1 the dye has propagated a unit length along each stream, which is depicted 
by solid lines in the top panel. Since all four streams are disconnected so far, the dynamic 
tree still consists of four disconnected branches, each of which corresponds to a colored 
stream of unit length. At step 2 the streams a and b merge. This is reflected in the 
dynamic tree, where the nodes a and b are now connected into a single cluster. Notice 
that the leaves a and b are not directly connected in the static tree; this connection 
reflects a special property of the dye's downstream propagation. 

At step 3 stream c reaches stream a. Since stream a by that time is already merged 
with stream b, we say that the stream c merges with the cluster of a and b; this is 
reflected in the dynamic tree in the lower panel for this step. Hence, at step 3 there exist 
two connected clusters of the colored flux: one cluster is formed by the streams a, b, 
and c, while stream d alone forms the second cluster. Finally, at step 4, all the colored 
fluxes merge together. The conventional representation of both static and dynamic trees, 
which does not show the link lengths, is given in the two rightmost panels. 

This example shows that the dynamic tree Td can be very different from the corre- 
sponding static tree T5. We notice in particular that in this example the static tree is 
a tree with no side branching; it has the largest possible Horton-Strahler order, 1] = 3, 
for a tree with four leaves. At the same time, the dynamic tree exhibits exhaustive 
side-branching; accordingly, it has the smallest possible order, Q = 2, for a four-leaved 
tree. 

3.2 Realistic example 

Here we illustrate the dynamic tree for an order-3 subbasin of the Noyo basin; this basin 
is located in Mendocino County, California, USA, and is described by Sklar et al. [2006]. 
The stream network for this subbasin is shown in Fig. HI its fifteen sources are marked 
by numbers 1 to 15 and fourteen stream joints by letters a to n. The static tree Ts for 
this stream network is shown in Fig. [5^; it has the Horton-Strahler order f2 = 3. 

The time-oriented dynamic tree T/) is shown in Fig. [5]d against the time axis (on the 
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ordinate); notice that time can also be interpreted as the distance traveled by the dye 
from each source. This interpretation has a direct connection to the metric properties of 
the basin and we will use it in the subsequent analysis. The order of the djTiamic tree is 
Q = 4. The letter and number marks in Fig. [5] match those in Fig. HI 

Four snapshots of the dye propagation — at times t = 1, 20, 39, and 60 — are shown 
in Fig. [6l In this example, the dynamic tree shows a larger degree of side-branching 
compared to the static tree; this larger degree is reflected in its larger HS order. We shall 
see in other realistic examples, further below, that this seems to be the case for most 
actual river networks. 

4 Stream vs. hillslope networks 

In an actual landscape, channels are initiated when the area upstream suffices to cre- 
ate a sustainable source of streamfiow and this source imprints a permanent channel on 
the terrain. Although these channels are typically detectable by field observations, the 
extraction of the channel initiation points, or "channel heads," from Digital Elevation 
Models (DEMs) has been a subject of intense study. Most commonly, channels are as- 
sumed to be initiated when the upstream area, or area times a typical slope, exceed a 
given threshold; the parameters of such relationships are field-calibrated. More recently, 
the availability of high-resolution, 1-m elevation data from Light Detection and Ranging 
(LIDAR) instrumentation has initiated a new generation of methodologies for the auto- 
matic detection of channels as "edges" or "features" in the terrain [e.g., Lashermes et 
al, 2007; Passalaqua et al, 2009]. 

The channelized paths, i. e. the branches of the river network, are not the only parts 
of the basin by which water or other fluxes — e.g., sediments, nutrients, or pollutants — 
are transported downstream. The unchannelized part of the basin, often called zero-order 
basins or hillslopes, is drained by pathways that have their own topology. In this work, 
we extract (i) stream networks from DEMs by using a critical threshold area Ac, and (ii) 
hillslope networks by assuming that Ac is as small as the DEM resolution. 

Clearly, each stream network is a part of the corresponding hillslope network. For a 
generic river basin, though, the total length of channelized paths is much smaller than 
the total length of unchannelized paths. In the present study, we assume that stream 
networks reflect the properties of channelized paths, while hillslope networks reflect the 
properties of unchannelized paths. The study of unchannelized-path topology below will 
show that it is quite different from the topology of the channelized paths. 
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Construction of stream and hillslope static trees is illustrated in Fig. [71 Figure [7^ 
shows a small part of a river basin; it is divided into 16 square regions called pixels. The 
well-defined streams occupy some of the pixels (shaded squares), the rest of the pixels 
(white squares) represent hillslopes, i.e. unchannelized parts of the basin. 

The elevation data can be used to figure out the flux direction from each pixel, whether 
stream or hillslope; this direction is depicted by arrows in Fig. [7)d. We assume that there 
is a unique flux direction away from each pixel; at the same time, fluxes can reach a 
given pixel from more than one other pixel. This property allows one to represent the 
directional information by a tree, which is shown in Fig. [7b. Solid nodes and solid lines in 
the figure represent the stream pixels and the stream flow respectively, while open nodes 
and dashed lines represent the hillslope nodes and hillslope flow. 

The final step in creating the static tree of this subbasin is to remove the linear 
segments (chains), that is to remove the nodes with only two connections (except the 
tree root). The resulting static hillslope tree is shown in panel (d). The static stream 
tree is obtained from the hillslope tree by removing the dashed links that represent 
unchannelized paths, and removing the remaining chains; the stream tree for our example 
is shown in panel (e). 

5 The dynamics of hierarchical aggregation 

The consecutive merging of river streams discussed in the previous section is a special 
case of a general phenomenon of hierarchical aggregation. This phenomenon is also called 
inverse cascading, and it can be described as follows. 

Consider a process that starts at time t = with N individual particles, which can be 
considered as clusters of unit mass. As time evolves, the clusters start to merge with one 
another, according to a set of suitable rules, thus forming consecutively larger clusters. 
We assume that only two clusters can merge at the same time; thus after each merging 
the number of clusters decreases by one. The process continues until all particles have 
been merged into a single cluster of mass N. The evolution of the above process can be 
described by a time-oriented binary tree, whose leaves correspond to the initial particles, 
the root to the final cluster of particles, and each internal node to the merging of a 
particular pair of clusters. 
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5.1 Examples 

Among the many instances of the above general aggregation scheme, we mention here 
the following four. 

Percolation. In the site percolation process on an L x L lattice, the initial N — L"^ 
particles correspond to the sites of the lattice, while clusters correspond to connected 
patches of occupied sites that are formed during the percolation process [Zaliapin et 
al., 2005]. In fact, the same scheme can be apphed to bond percolation, as well as to 
percolation on grids in higher dimensions. 

Billicirds. Elastic billiard on a rectangular table can be used to model gas dynamics in 
two dimensions (2-D). Here the initial particles are the N billiard balls (gas molecules) at 
time t — Q. Each of the balls is assigned an initial position and velocity. The clusters at 
time A are formed by balls that have collided during the time interval [0, A] [Gabrielov 
et al, 2008]. Formally, two balls are called A- neighbors if they coUided during the time 
interval [0, A]. Each connected component of this neighbor relation is called a A- cluster. 
Notice that within an arbitrary A-cluster each ball has collided with at least one other 
ball during the time interval [0, A]. In other words, a A-cluster is a group of balls that 
have affected each other's dynamics during the time interval of duration A. The mass of 
each cluster is simply the total number of balls within that cluster. Upon many collisions 
of the balls, the whole system will be composed of clusters of different sizes. As time 
evolves, the number of clusters will decrease and their mass increase. 

The same scheme can be applied to a system of particles that interact according to 
some potential f/(x). Bogolyubov [1960] suggested that when the interaction of particles 
is short-ranged, the system can be decomposed into finite clusters so that during some 
random interval of time, each cluster moves independently of other clusters as a finite- 
dimensional dynamical system. After this time interval, the system can be decomposed 
again into other dynamically independent clusters and so on. This type of dynamics 
is called cluster dynamics and Sinai [1974] showed analytically that it exists in a one- 
dimensional (1-D) system of statistical mechanics. Numerical results of Gabrielov et al. 
[2008] describe the presence and various properties of cluster dynamics in a 2-D system 
of hard balls. 

Phylogenetic trees. Probably the best-known application of hierarchical aggregation 
is in constructing phylogenetic trees that describe the evolutionary relationships among 
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biological species [Maker, 2002]. Here, a node corresponds to a set of species. Two species 
are connected if they have a direct common ancestor; the link length from a species to 
its direct ancestor equals the time it took to develop the descendant species from that 
ancestor. 



River transport. The example of interest to us here is the downstream transport 
along a river network. In this case, the initial particles are the environmental fluxes at 
the sources of the network, and clusters are formed by consecutive merging of the streams 
down the river path. That is, new clusters are formed when fluxes from upstream merge 
at the stream junctions. This scheme of describing dynamics along a static tree was 
considered in detail in Section [Sj albeit without referring to hierarchical aggregation. 



5.2 General set-up 

Hierarchical aggregation can be described in great generality by using the framework of 
nearest-neighbor clustering in a metric space. Specifically, consider a set S with distance 
d{a, h) for a, 6 G S; the elements of the set will be called points. The distance d{A, B) 
between two subsets of points A = {cii}i=i,...,Ar^ and B = {6i}i=i,...,ArB from § is defined 
as the shortest distance between the elements of the sets: 

d(A,B)= min diai^hj). 

l<i<NA,l<j<NB 

Nearest-neighbor clustering is a process that combines points from § into consecutively 
larger subsets, called clusters, by connecting at each step the two nearest clusters; this 
process can be described by the nearest-neighbor spanning tree T. Specifically, consider 
N points c° G §, i = 1,...,N with pairwise distances = (i(c°, c°). These points, 
considered as clusters of unit mass [rrii = 1), form leaves of the time-oriented tree T. 
The first internal tree node is formed at the time ti = min^j d^j by merging two closest 
points c° and c°, with [i*,]*) = argmin (i° • , where argmiUj^ /(i, j) is defined as a pair 
such that f{i*,j*) = minjj f{i,j). This merging creates a new cluster of two 
points, with a mass of rrii -\- rnj = 2. Hence, at time ti, there exist — 1 clusters: N — 2 
clusters with unit mass and one cluster of mass m = 2. 

We can now reindex the clusters so as to work with clusters cj, i = 1, . . . , N — 1; 
their total mass is Xlila^ nii = N and pairwise distances are djj = d{c}, c]). The second 
internal node of tree T is formed at time t2 = min^j djj > ti by merging the two closest 
clusters from the set {Cj^}j=i^..._iv-i- Thus, at time t2 we have N — 2 clusters such 
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that their total mass is and pairwise distances are dfj = d{cf,Cj). We continue in 
the same fashion, so the k-th internal cluster, for l<k<N — 2, is formed at time 
tk = minjj d'^j > tk-i, and at that time we have (A^ — k) clusters c'^ , i = 1, . . . , N — k with 
masses rrii such that X^^ITl'^ — Finally, at time t^^i we create a single cluster of 
mass that combines all points c^; this cluster forms the root of the tree T. 

Consider two nodes a and b from the nearest-neighbor tree and let and be their 
time marks; recall that the tree is time-oriented by the definition of the successive times 
tk = minjj c/^j- > tfc-i at which the cluster mergers occur. The ancestors of a node are its 
parent, the parent of that parent, and so on, all the way to the root. Clearly, the time 
mark for an ancestor is larger than that of a descendant. The nearest common ancestor 
p of nodes a and b is their common ancestor with the minimal time mark tp. 

The distance u{a, b) along the the nearest-neighbor tree is defined as the maximum 
of the values u{a,p) = tp — ta and u{b,p) = tp — tf,. This distance satisfies two of the 
usual distance axioms, symmetry and strict positivity, but the triangle inequality can be 
replaced by a more stringent one, namely 

u{a, b) < max [u{a, c), m(c, b)], 

which holds for any three nodes a, b and c. Such a distance function is called an ultra- 
metric [Rammal et ai, 1986; Schikhof, 2007]. Ultrametric spaces have many peculiar 
properties; for instance, one can rename any triplet a, b, c of nodes in such a way that 

u{a, c) = u{b, c). 

These unusual properties give ultrametric spaces considerable flexibility in applications, 
and point sets connected via nearest-neighbor clustering are a representative example of 
such spaces. 

In the billiard example of Section 15. H the space § is the set of billiard balls and 
the ultrametric distance function u{a, b) equals the time before the first collision of the 
balls a and b. Naturally, their distance depends on the initial positions and velocities of 
the two balls a and b, but it is affected by the global billiard dynamics: our two balls 
may be set to collide at a given time t* in the absence of other balls, but may be hit by 
some other ball at time t < t*, thus postponing the collision of a with b. 

In our river transport problem, the space § is the set of all river sources. The ul- 
trametric distance u{a, b) between two sources is defined as the time necessary for the 
corresponding fluxes injected into these two sources to meet down the river path. If the 
static river geometry is described by the tree Ts — and we assume, as previously stated. 
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that fluxes move with unit speed downstream — the traditional distance d{a, h) between 
two sources equals the maximal length along the tree to their nearest common parent in 
Ts- The nearest-neighbor spanning tree of hierarchical-aggregation theory becomes what 
we called so far, in the context of river transport, the dynamic tree Td- As previously 
stated, this dynamic tree differs, in general, from the static tree Ts and depends not only 
on the topology of the latter, but also on the actual length of the links. If the velocities 
vary in time or space, then the spanning tree Td will depend on the specific dynamics of 
the processes operating on the static tree. 

To better understand transport on river networks, we will elucidate in the next section 
the connection between the statistical properties of Ts and those of T^. 

6 Analysis of drainage networks 

In this section we quantify similarities and differences between the branching topology 
of static and dynamic trees, both at the stream and hillslope network scale. 

6.1 Data description 

We have analyzed three river basins: Upper Noyo (Mendocino County, California, USA), 
Tirso (Sardinia, Italy), and Grigno (Trento, Italy). Infomration about the physiographic 
and geologic characteristics of these basins can be found in, respectively, Sklar et al. 
[2006], Pmna et al. [2004], and Guzzetti et al. [2005]. The available DEMs were at a 
resolution of 10x10 m^ for the Noyo basin, 30x30 m^ for the Grigno basin, and 100x100 
m^ for the Tirso basin. Since the focus of this study is not the extraction of the most 
accurate river network from the available DEMs, we felt comfortable adopting a simple 
criterion for channel initiation as 100 pixels for all basins. Our main conclusions about 
the comparison between the static and dynamic trees would not be affected by changing 
the critical threshold within reasonable ranges. 

The static trees we extracted from these DEMs for the three stream networks are 
shown in Fig. [HI Using the procedure described earlier, we also extracted the static trees 
for the hillslope networks, which drain every pixel of a basin, by using a steepest gradient 
algorithm. The corresonding dynamic stream and hillslope trees were then constructed 
for each basin, assuming a constant unit speed of downstream propagation for the fluxes. 
Thus, we analyzed four different kinds of tree — static stream, dynamic stream, static 
hillslope, and dynamic hillslope - for each basin. 
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6.2 Self-similar properties 

Figures [91 and [TOl show the distributions of the number A^^, average magnitude Mr, and 
the average number Cr of hnks for branches of order r. The results in Fig. [9] refer to the 
stream trees; the results in Fig. [10] to the hillslope trees. 

Despite the small-sample fluctuations, the figures demonstrate a large degree of con- 
sistency among the branching indices for the trees from different classes. All considered 
branching statistics are closely approximated by the Horton laws. Moreover, the results 
suggest that the relationship holds in all the considered cases. Furthermore, we ob- 
serve that the values of the stream ratios for static trees are higher than the corresponding 
values for dynamic trees; and the values of the stream ratios for stream trees are smaller 
than the corresponding values for hillslope trees. 

The only indices that considerably deviate from the Horton laws at higher orders 
are Cr (average number of nodes within an order-r branch) for the Noyo basin and this 
warrants special investigation in the future. Apart from this discrepancy, overall we 
conclude that the four classes of trees, dynamic vs. static and stream vs. hillslope, can 
be closely approximated by the Tokunaga SSTs. 

6.3 Phase transition in hierarchical dynamics 

Here we ask the question as to whether the river network connectivity (in terms of ele- 
ments of the network participating in transport) exibits a phase transition akin to those 
found in other systems. Figure [TT] shows the fractional magnitudes rrii/N of the branches 
in the dynamic trees (stream and hillslope trees of the three river basins) as a function 
of the distance d traveled by the dye. Recall that this distance can also be interpreted as 
the time t when the node was created by merging of upstream branches. Altogether we 
consider six cases; in all of them one observes the following scenario. We start at distance 
d = (or time t = 0) with A^ branches (clusters) of unit magnitude corresponding to 
the most outer nodes of the transport tree. As distance increases (time evolves), the 
number of clusters decreases while their magnitudes become larger and exhibit promi- 
nent variability. In particular, at the small distances the maximal magnitude increases 
exponentially with distance; this is reflected by an approximately linear form of an upper 
envelop of the points in the figures (the envelop is not shown). Furthermore, we notice 
that at the small distances (times) the magnitude distribution is "continuous" in a sense 
that it does not have significant gaps. However, at some critical time t* (translated here 
to distance d* for easier interpretation), the distribution undergoes a serious qualitative 
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change: a prominent maximal cluster appears, such that its magnitude becomes signifi- 
cantly larger than that of the second larger cluster. Moreover, while the magnitude of the 
largest cluster keeps growing, the rest of the distribution is fading off so after some time 
all clusters present ai d = merge with the largest cluster. An interesting observation 
is that at the critical distance d* the magnitude of the largest cluster is just about 10% 
of the total magnitude of the system. Notably, this number is universal for all the 
considered examples. 

Figure [I2] shows the magnitude distribution of the clusters that existed when the dye 
traveled a given distance d. The analysis is done for the critical distance d* and a smaller 
distance d ~ d*/2; they are both indicated by vertical lines in Fig. [TTl In all six cases, we 
see that the magnitude distribution at the smaller distance (squares) has an exponential 
tail, while at the critical distance (circles) it is a power law. Recall that, in a log-log plot, 
power-law behavior shows up as a straight line, while exponential behavior becomes a 
convex curve. This change indicates that a phase transition occurs at the distance d*. 

This phase transition is further illustrated in Fig. [131 which shows three snapshots of 
the dye propagating down the Noyo basin. The distances traveled by the dye at these 
snapshots are marked by vertical lines in Fig. [HI the figure shows the number of clusters 
(dotted line) and the magnitude of the largest cluster for the Noyo dynamic tree (solid 
line), as a function of downstream propagation distance. 

The values of the six critical distances d* shown in Fig. [TT] vary over two orders of 
magnitude and depend strongly on the particular network being analyzed. Nevertheless, 
we notice a very good power-law fit for the value of d* in terms of the average link length 
L of the corresponding static tree (see Fig. [T5l) : 

d*^3.5L. (10) 

This relationship can be interpreted as follows in terms of the transport on river networks: 
the giant cluster of connected streams is formed when each flux traveled approximately 
3.5 links downstream from a source. We conjecture that: (a) this is a universal property 
of downstream transport on Tokunaga trees with rich branching, i.e. Tokunaga SSTs 
with c > 1 in Eq. ([5]); (b) the coefficient of proportionality in (ITOl) may depend on the 
Tokunaga parameters, but only weakly; and (c) this coefficient is larger than or equal to 
2 for any binary tree. An in-depth investigation of this issue is left for future study. 
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7 Concluding remarks 



7.1 Summary 

This study focused on the statistical description of environmental transport on self-similar 
river networks. We approached the problem by considering downstream transport on such 
a network as a particular case of nearest-neighbor hierarchical aggregation; the so-called 
ultrametric induced by the branching structure of the river network provides the distance 
function with respect to which the downstream flow gives rise to clusters that decrease 
in number and increase in size with time. 

We described the static topological structure of a drainage network by the type of 
tree structure that goes back to the pioneering studies of Horton (1945), Strahler (1957) 
and Shreve (1966), and referred to it as a static tree, to distinguish it from the associated 
dynamic tree. This novel concept introduced herein describes downstream transport 
along the static tree. 

We studied the statistical properties of both static and dynamic trees using the 
Horton- Strahler (HS) and Tokunaga branching taxonomies [see Horton, 1945; Strahler, 
1957; Tokunaga, 1978]. Using three river networks — the Noyo, Grigno and Tirso — 
we showed that both static and dynamics trees can be well approximated by Tokunaga 
self-similar trees (SSTs). The HS and Tokunaga parameters of these two types of trees 
differ significantly, though, for each of the three basins. This difference supports the rel- 
evance of the dynamic tree concept; its parameter values depict important properties of 
the transport on a given river network that are not captured by the conventional, static 
tree. 

A striking result of this study is the phase transition we found in river network 
dynamics: as one fills an empty river network through its sources, or injects a dye into a 
water-filled one, the number of clusters of connected nodes decreases and the size of the 
largest cluster increases, until a dominant cluster of connected streams forms. During 
this process, the time-dependent size distribution of the connected clusters changes from 
an exponential to a power-law function as the critical time approaches. 

This phenomenon, which may seem rather unexpected in the present, hydrological 
setting, can be better understood within the framework of complex networks. This 
framework has been explored in many natural and socio-economic settings, ranging from 
the functioning of a cell to the organization of the internet [Albert and Barahasi, 2002]. 

The mathematical theory of complex networks considers a group of nodes that can 
be connected with each other according to some problem-specific rules, thus forming a 
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graph. In the simplest case, the node connections are independent of each other and can 
be specified by the probabihty p that two randomly chosen nodes are connected. There 
exists a critical value pc such that for p < Pc the network consists of isolated clusters, 
while a single giant cluster appears as p crosses pc, and spans the entire network. The 
appearance of this giant cluster is remarkably reminiscent of infinite-cluster formation in 
percolation theory [Stauffer and Aharony, 1994]. Albert and Barabasi [2002] provide a 
comprehensive review of parallels and differences between complex-network theory and 
percolation theory. 

It readily follows from the analysis of Section [3] that the transport on river basins 
fits rather naturally the complex-network paradigm. Formally, each river source is repre- 
sented by a node and two streams are considered to be connected when their respective 
fiuxes join downstream. This is exactly the scheme we used to define a dynamic tree, 
with the only difference that we ignored the node connections within already formed 
clusters. This difference does not affect the process of cluster formation, so all the results 
of complex-network theory do apply to the envirodynamics of river basins. 

There is an important difference, though, between complex networks in general and 
the dynamic trees considered in this study. Our dynamic trees, unlike general networks, 
are time-oriented, i.e., their nodes can be ordered according in "time" or with respect 
to a "distance" parameter. The ultrametric distance along such trees satisfies a stronger 
triangle inequality than ordinary distance (see Section 15. 2p . Spaces equipped with an 
ultrametric u, instead of a traditional distance d, have therefore interesting properties [e.g. 
Schikhof, 2007]. As shown in Section O hierarchical aggregation via nearest-neighbor 
clustering provides a common framework for many apparently different processes — e.g., 
billiards, river transport, and percolation — in the setting of ultrametric trees, and thus 
may provide novel insights into these processes. 

In percolation models, the cluster-size distribution at phase transition is given by a 
power law, whose index is a function of the system's dimension alone. In our three river 
networks, this index differs from the one to the other, and from the river to the hillslope 
network for the same basin. In our hierarchical aggregation on dynamic trees, different 
clustering rules may correspond to different effective "dimensions" of the system. At 
the same time, it is known that the critical percolation indices are universal for systems 
in high dimensions [Hara and Slade, 1990] and trees are a simple model for infinite- 
dimensional systems [Albert and Barabasi, 2002]. Thus, one expects to see the same values 
of the critical indices when working with percolation on a tree. From this perspective, 
the fact that our critical exponents vary from basin to basin, and from river to hillslope 
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trees, still needs to be understood. 

7.2 Discussion and further work 

In this study we considered only the simplest clustering rules for the river streams: two 
streams belong to the same cluster if there is a connected path from one stream to another 
along the river network. This approach is patterned after percolation studies and allows 
for a straightforward treatment. It might however result in a situation when two streams 
belong to the same cluster despite the fact that the respective fluxes are not mixed yet 
(think of two short streams that merge with a spatially extended cluster at about the 
same time). Formulating a physically more appropriate set of clustering rules might yield 
more realistic results for a wealth of river networks with differing properties. 

So far, we only investigated dynamic trees that have the same set of leaves as the 
corresponding static tree; this corresponds to injecting a flux through the sources. At the 
same time, it might happen that a flux of interest is injected into an internal node, e. g., 
an industrial pollutant from a plant or nutrient production from a local biotic activity. 
Such situations can be easily modeled by considering a dynamic tree whose set of leaves 
samples the entire river network. 

To construct a richer theoretical framework for transport on river networks one may 
also model the transport along real and synthetic networks by using Boolean delay equa- 
tions (BDEs). In BDEs, the discrete state variables describe the flux through the river 
branches; naturally, the rules for updating these variables inherit the child-parent rela- 
tionship of the stream's static tree. The parent variables are updated based on the values 
of the children variables, after delays that correspond to time of flux propagation from a 
child to its parent. Ghil et al. [2008] reviewed recently BDEs and their applications to 
climate and earthquake modeling. We expect such modeling to shed further light on the 
complex and important problems of river transport. 
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(a) River network (b) Corresponding tree 



Figure 1: Tree representation of a river network: (a) hypothetical river network; and (b) 
its representation by a binary tree. The network sources and the respective tree leaves 
are marked by the same letters in both panels. The figure also illustrates the terminology 
used in our river transport study. 
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(a) Horton-Strahler orders (b) Tokunaga indices 

Figure 2: Example of (a) Horton-Strahler ordering, and (b) Tokunaga indexing. 
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Figure 3: Constructing a dynamic tree. The initial static tree and the final dynamic tree 
are shown in the rightmost pair of panels. The dynamic tree refiects the propagation of 
a fiux from leaves to the root of the static tree at a constant velocity. The top row of 
panels shows the static tree at different steps of this process; for visual convenience we 
explicitly show the static tree's link lengths. The bottom row shows the corresponding 
phases of the dynamic tree. The top leftmost panel indicates the lengths of the links in 
the static tree; each step in the figure takes one time unit, that is the fiux propagates 
one unit of length downstream. See Section 13.11 for details. 
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Figure 4: Stream network for an order-3 subbasin of the Noyo river basin, Mendocino 
county, California. Sources are marked by letters, stream merging points by numbers. 
The same marks are used in Figs. [5] and [6] that show the static and dynamic trees for 
this stream. 



(a) (b) 




n 

Figure 5: Static and dynamic trees for the Noyo subbasin of Fig. HI (a) Static tree Ts 
and (b) dynamic tree Td. Letter and number marks are the same as in Fig. HI 
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(a)t=1 (b)t^20 (c)t^39 (d) t=60 




Figure 6: Three snapshots of the evolution of the dynamic tree (heavy sohd hues) on the 
static tree (hght sohd hues) for the stream of Fig. HI Letter and number marks are the 
same as in Fig. |H 
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Pixellzed river basin 



(b) Flux directions 




(c) Hillslope tree with linear (d) Hllislope tree (0 = 3) (e) Stream tree {U =2) 
segments (chains) 



Figure 7: Construction of a static tree tliat represents tlie topology of liillslope (unclian- 
nelized) and stream (cliannelized) drainage patlis. (a) Pixelized river basin; the shaded 
pixels (cells) correspond to the stream location (solid line), the white pixels - to the 
valleys or hillslopes. (b) Flux direction obtained from the elevation data, (c) Tree that 
describes the drainage topology: solid nodes and links correspond to the stream pixels 
and stream flow; open nodes and dashed links - to the hillslope pixels and hillslope flow. 
Notice that this tree contains several purely hnear segments, with no branching, (d) The 
same tree, from which the linear segments have been removed: it describes the topology 
of both hillslope paths and stream paths, and is referred to as the hillslope tree, (e) The 
subset of the tree in panel (d) that describes the topology of the stream paths only; this 
tree is referred to as the stream tree. 
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(A) Noyo stream 




(B) Grigno stream 




8SS 
3fV 



(C) Tirso stream 



























10,000 meters 







Figure 8: Static trees for the stream networks of the three basins analyzed in this study, 
a) Noyo, Mendocino County, Cahfornia, USA; the outlet is marked by a ball; b) Grigno, 
Trento, Italy; c) Tirso, Sardinia, Italy. See Section 16.11 for details of channel initiation. 
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Figure 9: Branching statistics for the stream trees of Noyo, Grigno, and Tirso basins. 
Number and average magnitude Mr for static (panel a) and dynamic (panel c) trees 
and average number of links within a branch for static (panel b) and dynamic (panel 
d) trees. 
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Figure 10: Branching statistics for the hillslope trees of Noyo, Grigno, and Tirso basins. 
Number A^^ and average magnitude for static (panel a) and dynamic (panel c) trees 
and average number of links within a branch for static (panel b) and dynamic (panel 
d) trees. 
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Figure 11: Fractional branch magnitudes irii/N as a function of the distance di traveled 
by the dye at the branch creation instant. Notice that the distance di can be interpreted 
as the time ti necessary to create the branch, a) Noyo stream; b) Noyo hillslope; c) 
Grigno stream; d) Grigno hillslope; e) Tirso stream; f) Tirso hillslope. 
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Figure 12: Distribution of branch magnitudes at the critical distance d* (balls) and 
at an earlier time (squares) in dynamic trees for the three basins, a) Noyo stream; b) 
Noyo hillslope; c) Grigno stream; d) Grigno hillslope; e) Tirso stream; f) Tirso hillslope. 
Each panel shows two distributions, the corresponding distances are depicted by vertical 
lines in Fig. [TTl Notice that the value of the critical distance d* can be interpreted as the 
critical time t* necessary to create the critical cluster. The downward deviations from 
the pure power laws are due to the finite-size effect. 
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Figure 13: Transport down the Noyo stream network. Three snapshots of flux propa- 
gation from the stream sources to the outlet, at (a,d) d = 200, (b,e) d = 500, and (c,f) 
d = 1000. Panels (a)-(c) show the entire Noyo basin, while panels (d)-(f) zoom onto an 
order-4 subbasin located in the lower right part of the entire basin. See also Fig. [TH 
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Figure 14: Cluster evolution for the Noyo downstream flux transport: number (dotted 
line) and largest-cluster size (solid line). Vertical lines correspond to the three snapshots 
in Fig. I 
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Average link length, [m] 

Figure 15: Critical distance d* as a function of the average link length L for the six 
dynamic trees shown in Fig. [TTl The line in the figure corresponds to d* = 3.5 L. 
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